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ABSTRACT 

Gough & Mclntyre have suggested that the dynamics of the solar tachocline are 
dominated by the advection-diffusion balance between the differential rotation, a large- 
scale primordial field and baroclinicly driven meridional motions. This paper presents 
the first part of a study of the tachocline, in which a model of the rotation profile below 
the convection zone is constructed along the lines suggested by Gough & Mclntyre and 
solved numerically. In this first part, a reduced model of the tachocline is derived in 
which the effects of compressibility and energy transport on the system are neglected; 
the meridional motions are driven instead by Ekman-Hartmann pumping. Through 
this simplification, the interaction of the fluid flow and the magnetic field can be 
isolated and is studied through nonlinear numerical analysis for various field strengths 
and diffusivities. It is shown that there exists only a narrow range of magnetic field 
strengths for which the system can achieve a nearly uniform rotation. The results are 
discussed with respect to observations and to the limitations of this initial approach. 
A following paper combines the effects of realistic baroclinic driving and stratification 
with a model that follows closely the lines of work of Gough & Mclntyre. 



1 INTRODUCTION 



The solar tachocline is a thin shear layer located at the 
top of the radiative zone. It performs the dynamical tran- 
sition between the differentially rotating convection zone 
and the nearly uniformly rotating radiative zone (Brown et 
al., 1989). In a non-rotating fluid, the interface between a 
convective and a non-convective region is intrinsically com- 
plex, and supports dynamical phenomena such as overshoot 
and gravity-wave generation. Including a background ro- 
tation increases the complexity of the system by driving 
meridional motions at the interface (Mestel, 1953) and by 
changing the characteristics of convective motions and over- 
shooting plumes, both leading in particular to non-vanishing 
anisotropic stresses and differential rotation. Finally, the 
combination of strong shear and anisotropic small-scale mo- 
tion has been shown to lead sometimes to the enhancement 
of any seed magnetic field, and to intermittent (or even- 
tually quasi-periodic) magnetic phenomena which could be 
related to the observed 11-yr magnetic cycle. The position 
of the tachocline in this unique dynamical interface makes 
it one of the most complex, and as a result, one of the most 
interesting regions of the sun. 

In recent years, the tachocline has been one of the prin- 
cipal targets of helioseismic observations of the angular ro- 
tation profile of the sun. The determination of the thickness 
of the tachocline is essential to understanding its dynam- 
ics. Assuming a variation profile for the angular velocity 
across the tachocline, throughout the radiative interior and 
the convection zone, artificial frequency splittings can be re- 



constructed which are then fitted to the observations. The 
most recent estimate of the thickness of the tachocline A 
using this method is the following: A = (0.039 ± O.O13)R,0 
(Charbonneau et al.,1999). An independent method is pro- 
posed by Elliott & Gough (1999): they suggest that the 
discrepancy between the observed sound-speed profile and 
that obtained from Standard Solar Models could be related 
to a chemical composition anomaly caused by the mixing 
of helium within the tachocline layer. By including an ad- 
ditional mixing layer below the convection zone, and com- 
paring the predicted sound speed of this new model to the 
observations, they calibrate the thickness of the tachocline 
to be approximately 0.02 R©, assuming the tachocline to be 
strictly spherical. If the tachocline is thinner but aspherical, 
this estimate is a measure of the asphericity. Note that ob- 
servations of light-elements abundances at the surface also 
suggest that mixing below the convection zone is confined 
to a shallow layer (e.g. Brun, Turck-Chieze & Zahn, 1999). 

It has also been suggested that helioseismic observa- 
tions could be used to determine the latitudinal structure 
of the tachocline, in particular whether the width and posi- 
tion of the tachocline vary with latitude. Independent anal- 
yses by Gough & Kosovichev (1995) and Charbonneau et 
al. (1999) suggest that the base of the convection zone and 
the tachocline (respectively) may be prolate. Charbonneau 
et al. (1999) also report that no significant variation in the 
tachocline thickness can be deduced from the observations. 
Finally, the most recent observational results on the struc- 
ture of the tachocline concern its dynamical aspect: Howe 
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et al. (2000) suggest that a large-scale oscillation may be 
taking place across the tachocline with a period of about 
1.3 yr, although these results are contested (Basu & Antia, 
2001). If assumed to be related to Alfvenic torsional oscil- 
lations, this observation could be evidence for the presence 
of a magnetic field threading the tachocline with a radial 
component of about 500 G (Gough, 2000). 

Owing to the complexity of the dynamics of this region, 
only a handful of models of the tachocline have been pro- 
posed so far. Spiegel & Zahn (1992) pointed out that the 
shear across the tachocline must be associated with signifi- 
cant thermal fluctuations through a thermal-wind balance. 
These temperature fluctuations drive meridional motions 
and propagate the shear deeper and deeper into the radia- 
tive zone. In order to avoid this radiative spreading, Spiegel 
& Zahn proposed an anisotropic turbulence stress-model to 
suppress the latitudinal shear within a short vertical length- 
scale, the strong stratification of the radiative zone providing 
a natural explanation for the anisotropy of the small-scale 
flow. However, it was later argued that such a model can- 
not adequately describe the tachocline (Gough, 1997, Gough 
& Mclntyre, 1998). Broadly speaking, the strong restric- 
tion of turbulent motions to spherical shells by the back- 
ground stratification in the radiative zone leads to Reynolds 
stresses that transport not angular momentum but potential 
vorticity (see Garaud, 2001a, for instance). These Reynolds 
stresses would therefore drive the system away from rather 
than towards uniform rotation. Moreover, Spiegel & Zahn's 
model of the tachocline assumes that the tachocline is tur- 
bulent; it is not clear whether this is indeed the case. The 
tachocline is likely to be stable (Charbonneau, Dikpati & 
Gilman, 1999) or marginally stable (Garaud, 2001a) to hy- 
drodynamical linear shear instabilities. However, by assum- 
ing the existence of a relatively strong toroidal field, it was 
shown that MHD 2-D instabilities exist (Dikpati & Gilman, 
1999) and can have a significant influence on the redistribu- 
tion of angular momentum (Cally, 2001). 

One of the most natural models for the solar interior 
rotation, and in particular the quenching of the shear, is 
to assume the existence of a large-scale primordial field 
in the radiative zone. Indeed, as the magnetic diffusion 
timescale in the radiative zone is much larger than the rota- 
tion timescale, Ferraro's isorotation law (1937) ensures that 
the angular velocity is constant on field lines. Mestel & Weiss 
(1987) suggested that large-scale fields of amplitude as low 



as 10" 



10 G should be capable of suppressing most of 



the rotational shear deep in the solar radiative interior. Sim- 
ulations performed by Rudiger & Kitchatinov (1997) and 
MacGregor & Charbonneau (1999) looked at the interac- 
tion between a large-scale field with a fixed poloidal compo- 
nent and a purely azimuthal flow. These seemed to confirm 
this estimate, and indeed managed to reproduce the con- 
finement of the shear to a thin tachocline. However, as it 
was pointed out by Gough & Mclntyre (1998), the very low 
magnetic diffusivity of the radiative zone enables a strong 
nonlinear coupling between the thermally driven meridional 
flows described by Spiegel & Zahn and the large-scale in- 
terior field. This interaction is essential to the dynamics of 
the tachocline. Starting from this idea, Gough & Mclntyre 
(1998) presented the first model of the tachocline to take 
into account self-consistently the nonlinear interaction be- 
tween thermally driven meridional flows and a large-scale 



magnetic field. The complexity of the resulting model, how- 
ever, precluded the derivation of a complete solution. 

It is the purpose of this work to create a model of 
the tachocline which encompasses gradually more realis- 
tic physics, to reproduce eventually the idea proposed by 
Gough & Mclntyre. Only this type of gradual approach can 
provide an adequate procedure to isolate different dynami- 
cal phenomena and to study them individually before com- 
bining them into one complex system. In this first paper, 
the compressibility of the fluid is neglected (as well as the 
background stratification and the energy transport), in or- 
der to study specifically the nonlinear interaction between 
the magnetic field and the fluid motions. As a result, merid- 
ional motions must be driven artificially. This can be done 
through Ekman-Hartmann pumping on the boundary rep- 
resenting the interface with the convection zone, which has 
the main advantage of providing a well-known, contro lled 



flow pattern. This idea is discussed in detail in Section 6.3 



9 



Subsequent work will study the effects of stratification as 
well as advective and diffusive heat transport on the system, 
which can then provide a realistic quantitative description 
of a thermally driven flow and ultimately of the tachocline. 

The initial model studied in this first paper, is presented 
in Section ^: the assumptions and boundary conditions are 
discussed and applied to the system, and the numerical pro- 
cedure for the solution of the governing equations is outlined. 
In Section ^, the numerical procedure is tested by study- 
ing a non-magnetic case. The simulations are compared to 
the well-known Proudman-Stewartson solution for an in- 
compressible fluid between two concentric rotating spheres 
(Proudman, 1956, and Stewartson, 1966). The results of the 
simulations in the magnetic case are presented in Section 
It is shown that the effects of angular-momentum ad- 
vection and poloidal-field advection by the meridional flow 
are essential to the dynamics of the system, thereby em- 
phasizing the importance of correctly taking both into ac- 
count in realistic models of the tachocline. In particular it 
is shown that, contrary to common expectations, increas- 
ing the amplitude of the magnetic field does not necessarily 
reduce the amount of shear in the radiative zone. Instead, 
the simulations presented in Section ^ reveal the existence 
of three typical regimes, with low-, high- and intermediate- 
field-strength; only the intermediate regime is able to re- 
produce the qualitative features of the tachocline. Section 
fA presents two boundary-layer analyses which apply locally 
to different field configurations; the analytical solutions are 
compared to the numerical simulations. Finally, the simula- 
tions are analysed and discussed in Section n with respect 
to previous work and to observations. 



2 THE MODEL 

The principal aim of this study is to try to reproduce the 
observed rotation profile of the sun in the region below the 
convection zone, and in particular 

(i) the sharp transition (within a thin shear layer of thick- 
ness no larger than 4 % of the solar radius) between the lat- 
itudinal shear observed at the base of the convection zone 
and the uniform rotation of the radiative interior. 

(ii) the value of the interior angular velocity Q c (observed 
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to be about 93% of the surface equatorial angular velocity 

ficq ) ■ 

Moreover, the model should also aim to reproduce the 
confinement of the material mixing to the region of the 
tachocline, as suggested by helioseismic measurements (El- 
liott & Gough, 1999) and surface light-elements abundances 
(Brun, Turck-Chieze & Zahn, 1999). 



2.1 The assumptions 

The model presented in this paper essentially follows the 
idea proposed by Gough & Mclntyre (1998) by assuming 
the existence of a large-scale primordial field in the radia- 
tive zone and looking at its nonlinear interaction with fluid 
motions (both azimuthal and meridional) near the bottom of 
the convection zone in particular. In order to do so, nonlinear 
MHD equations are solved, in a first instance assuming that 
the fluid is incompressible (this paper) and in subsequent pa- 
pers in an anelastic approximation of a compressible fluid. 
This gradual increase in the complexity of the phenomena 
studied allows the independent study of the nonlinear in- 
teraction of the field and the meridional flow on the one 
hand, and the effects of heat transport and stratification 
on the other hand. The principal shortfall of the incom- 
pressibility assumption is that the mechanisms driving the 
meridional motions as described by Spiegel & Zahn (1992) 
are absent from this first model; instead, they are replaced 
by a judicious choice of boundary conditions which lead to 
Ekman-Hart man n pumping of a qualitatively similar flow 
(see Section |6.3[ ) . 

Two principal simplifications are performed on the sys- 
tem to reduce it to a numerically tractable problem: it is 
assumed that the system is axisymmetric and is in a steady 
state. Axisymmetry is a wholly unjustified, albeit standard 
assumption, which allows a drastic simplification of the 
MHD equations. The steady-state assumption, on the other 
hand, is reasonably well justified provided all dynamical 
timescales (for the rotation and the meridional circulation) 
are much shorter than the stellar-evolution timescale (which 
ensures that the background stratification does not vary dur- 
ing the fluid motion), the magnetic-braking timescale (which 
ensures that the external torques exerted on the system are 
small), or the magnetic-diffusion timescale (which ensures 
that the amplitude of the interior field does not vary signif- 
icantly during the fluid turnover time). 

In the following, it will always be assumed that the 
tachocline lies completely in the stably stratified radiative 
zone, which agrees reasonably well with observations (Char- 
bonneau et al. 1999). It will also be assumed that the ro- 
tation profile in the convection zone is independent of the 
tachocline dynamics so that the convection zone can be 
taken as a boundary condition of the model. Finally, it will 
be assumed that the region studied is located sufficiently 
far below the region of influence of the solar dynamo (Ga- 
raud, 1999) to avoid complications linked with the unsteady 
character of the dynamo field. As a result of this assump- 
tion, the region studied is also likely to be located below the 
overshoot region. 

Within this framework, the system studied is limited to 
the radiative zone, in a region located between two concen- 
tric spherical boundaries, the outer one near the base of the 



convection zone r ou t = r c , where r c = 0.7R©, the inner one 
at a radius ri n = 0.35r c . The inner core (occupying r < n n ) 
is cut out of the region of the simulation to avoid singulari- 
ties at the origin; the position of the inner boundary should 
not have a significant influence on the results. The numerical 
value of ri n is chosen in such way as to facilitate the compar- 
ison of these simulations with the work of Dormy, Cardin & 
Jault (1998), who solved a similar system of equations and 
boundary conditions for geophysical applications. 



2.2 The equations 

If the fluid motions are assumed to be axisymmetric and 
incompressible, the MHD equations in a frame rotating with 
angular velocity fi c = fi c e z are the following: 

2pfl c x u + Vp + pg h — j x B — pvV 2 u = , 

Vx(uxB)+?)V 2 B = 0, 

V-w = , 

V B = . (1) 

where the nonlinear terms in the meridional circulation 
(u ■ V)it are neglected (it is verified a posteriori that this is 
a good approximation - see Section |3,3| ), and g h — gh(r)e r 
is the gravitational acceleration associated with the hydro- 
static background. The angular velocity fi c = 2.84 x 10~ 6 
s _1 is chosen to take the observed interior angular veloc- 
ity of the sun; B is the magnetic field (with components 
(B r ,Bg,B,f,)), j is the electric current, and u is the circu- 
lation which has components (u r ,ue,u^). In order to derive 
these equations, the density p is assumed to be constant 
(po = lg cm -3 ) everywhere, as is the kinematic viscosity v 
and the magnetic diffusivity rj. This approximation is made 
to simplify the numerical procedure slightly, but can eas- 
ily be removed. Since, as it will be shown later, the typical 
values of the viscosity and magnetic diffusivity used in the 
simulations are several orders of magnitude larger than in 
the sun, it seems pointless to try and represent accurately 
the variation of these quantities in this first analysis. 
Using the following new system of units: 



[t] = 1/fic , [B] = Bo , [u] = r c O c 



(2) 



where Bo is the typical strength of the radial field in the 
interior, the equations become: 



2e z x u 
Vx(u x B) 



- Vp - g b + Aj x B + E V V u 



with 



A = 



poQ'ir'i 



E u = 



and ^ = fi^ 



(3) 



(4) 



where all the quantities are now dimensionless and e z is the 
unit vector parallel to the rotation axis. The Ekman num- 
bers E v and E v represent the ratios of the rotation timescale 
to the diffusive timescales, and the Elsasser number A is the 
ratio of the typical amplitude of the Lorentz force to that of 
the Coriolis force. The system is solved for (u r ,ue,u$) and 
(-B r , Be, Btj,) by solving the azimuthal component of the mo- 
mentum equation, the azimuthal component of the vortic- 
ity equation, the integrated induction equation and the az- 
imuthal component of the induction equation, together with 



© 0000 RAS, MNRAS 000, 000-000 



4 P. Garaud 



e = o 
ue=o 

dO/d6 = 
B e = 
B^ = 



r = r 

out c 
U r = Uq = 

a = a cz 



r in = 0.35 i 

U r = Uq = 

a = a k 

V 2 B = 



e-71/2 

u e =0 
d£2/d6 = 
B r =0 
B^ = 

Figure 1. Boundary conditions. 



the mass conservation equation and the solenoidal condition; 
those equations are 



2(e z x u)^ 
2[Vx(e z xu)] 
u r BQ — uqBt 
[Vx(ttxB)]^ 
Vw 
V B 



A(j x B) 4 , + E v {V 2 u)< t , , 
A[Vx(j x B)k + £„(V a u) 

-£„(V 2 B) , 
, 
o , 



(5) 



where u> = Vxu is the vorticity. 



2.3 The boundary conditions 

In this first analysis, the boundary conditions chosen for the 
system are the simplest possible ones that still guarantee the 
existence of a solution; as a result, they are not necessarily 
the most accurate representation of the dynamics of the sun 
(and in particular of the interface of the tachocline with 
the convection zone). The effects of the boundary conditions 
on the syst em, and possible improvements, are discussed in 
Section 



6.3 



The boundary conditions are summarized in Fig. |l|; the 
MHD equations presented in equations (^) are solved in a re- 
gion located between two spherical boundaries, impermeable 
and no-slip boundaries so that the radial and the latitudi- 
nal components of the circulation vanish on the boundaries, 
and the azimuthal component of the circulation is given by 
the rotation of the boundaries. As a result, the angular- 
velocity perturbation is Cl — Qi n — Q c on the inner boundary 
where Qi n is an eigenvalue of the problem (see below) and 
Cl — Q C z (8) — Q c on the outer boundary, where 



fi cz (#) = fieq(l — 0,2 COS 2 9 — 0,4 COS 4 6) 



(6) 



and fi eq = 1.07O, 



ai = 0,4, = 0.15 (according to the obser- 
vations presented by Schou et al. (1998)). On the equatorial 
plane, symmetry arguments determine the behaviour of the 
solutions: ug = 0, B r = 0, dCl/d9 — and B^ = 0. On 
the poles, regularity conditions impose ug — 0, Bg — 0, 



dCl/d9 = and B^ = 0. The bounding spheres are as- 
sumed to be imperfectly conducting and the region out- 
side [ri n ,r out ] supports no fluid motion, which implies that 
the field in that region satisfies V 2 B = 0, with conditions 
that B —* at infinity, and that the field structure be- 
comes purely dipolar as r — > 0. The amplitude of the ra- 
dial component of the field near the poles is fixed such 
that B r (ri n ,0 = 0) = Bq. This uniquely determines the 
boundary conditions for the magnetic field, as a relation be- 
tween B r and Bg on the one hand, and between B$ and 
dB^/dr on the other hand. As a result of these condi- 
tions, when the system is not rotating, the solution for the 
magnetic field is a purely poloidal dipolar structure with a 
field strength varying as 1/r 3 . When the system is rotating, 
these boundary conditions drive a meridional flow through 
Ekman-Hartmann pumping. This flow is used to mimic the 
baroclinicly driven flow predicted to occur in the tachocline 
(Spiegel & Zahn, 1992, Gough & Mclntyre, 1998), in order 
to study its interaction with the imposed magnetic field. 
The principal caveat of this approach is that the typical ve- 
locities of the meridional flow scales with E v and E n as an 
Ekman-Hartmann flow rather than as a thermally driven 
flow (i.e. with ug ~ r c fi and u r ~ <5eh^ where 5eh is the 
Ekman-Hartmann boundary- layer thickness). This discrep- 
ancy should influence quantitative predictions of the model 
only whilst leaving qualitative results unaffected. 

One of the main aims of this study is to be able to 
predict the angular velocity of the interior. In order to do 
so, the angular velocity of the inner core Sli n is treated as 
an eigenvalue of the problem, and an additional boundary 
condition is imposed on the system accordingly, namely that 
no net torque is applied to the region interior to r = n n 
(this is a necessary condition to guarantee a steady state). 
This condition, which determines uniquely the value of Qi n , 
is equivalent to requiring that the integral of the angular- 
momentum flux through the boundary vanishes: 



ir/2 



on 

pvr sin 9— h r sm9B r Bj, ) sinfldfl = 

or 



(7) 



One can notice through this expression the main justifica- 
tion for choosing conducting boundary conditions: when the 
region below n n is insulating, the radial component of the 
current (and thereby the toroidal field) must vanish on the 
boundary; in that case the angular-momentum flux through 
the boundary would be purely viscous. Since viscous stresses 
are normally thought to be negligible in the sun, such a 
model would be a poor representation of the dynamical 
structure of the radiative zone. 



2.4 Numerical method 

The system of equations ^ and the boundary conditions 
presented in Section [2.3| constitute a well-posed partial dif- 
ferential equations system with one eigenvalue. The numer- 
ical method chosen for the solution of this system is the 
following: 

(i) to write the system of equations with respect to a 
spherical polar coordinate system, respecting the symme- 
try of the problem as well as the regularity conditions on 
the poles, 

(ii) to expand the latitudinal dependence of the equations 
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into Fourier modes of 9, or equivalently Chebishev polyno- 
mials T n ( cos6), 

(iii) to solve the resulting ODEs and algebraic equations 
with respect to the independent variable r using a Newton- 
Raphson relaxation method. 

A more detailed description of the numerical procedure can 
be found in the work by Garaud (2001b). 

Manual mesh-stretching is preferred over automatic 
mesh-point allocation, as the latter was found to have 
low performance for boundary-layer thicknesses order of 
10 _6 (r ou t — n n ) or less. This low performance is probably due 
to the fact that as many Fourier modes are used, the allo- 
cation algorithms that were tried cannot choose adequately 
according the most relevant function against which the mesh 
should be stretched. When stretching the mesh manually 
(according to the width and positions of the boundary layers 
determined in Section |^), a minimum of a hundred points 
is allocated to each of the boundary layers. Typically, an- 
other 200-400 points are allocated to the interval outside 
the boundary layers, depending on the simulations. 



3 RESULTS IN THE NON-MAGNETIC CASE 

As a first step towards the resolution of the model presented 
above, a simplified system is solved in which the magnetic 
field is ignored. The hydrodynamics of a fluid between two 
concentric impermeable rotating spheres is a relatively well- 
studied problem, both analytically (since Proudman, 1956) 
and numerically (see for instance Dormy, Cardin & Jault, 
1998); these previous studies can be compared with the re- 
sults of the non-magnetic numerical solutions obt ained with 
the numerical procedure presented in Section 
accuracy and performance. 



2.4 to test its 



3.1 Numerical results 

As a example, the solution to the non-magnetic subset of 
equations (ft) with an Ekman number £„ = 8x 10 -6 is pre- 
sented in Fig. P When no magnetic field is present, for low 
enough Ekman number, fluid motion is dominated by Cori- 
olis forces everywhere except in two boundary layers near 
the spherical boundaries, and in a shear layer at the tan- 
gent cylinder. In the bulk of the fluid, angular velocity is 
more-or-less constant on cylinders: indeed, when viscosity is 
negligible, the fluid dynamics equations for the incompress- 
ible fluid reduce to 



(Sl c x u)$ = 



(8) 



which implies that u must be parallel to the rotation axis, 
and 



(Vx(Sl c x u)), =0 



(9) 



which implies that the angular velocity must be independent 
of 2 where z is the cylindrical coordinate that runs parallel to 
the rotation axis (Proudman, 1916, Taylor, 1921). Viscous 
effects are necessary in the boundary layers to ensure the 
smooth transition between the rotation profile in the bulk 
of the fluid and that imposed at the boundaries. 
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Figure 3. Analytical prediction for the interior angular velocity 
as a function of gap width 5, and results of simulations for (5 = 0.5 
and 8 = 0.65 for values of the Ekman number as shown. 



3.2 Comparison with analytical asymptotic 
analysis 

This structure was first studied by Proudman (1956) and 
Stewartson (1966), in the case where E v is asymptotically 
small. It is possible to show (see the Appendix) that in this 
limit the interior angular velocity is uniquely determined 
by the size of the gap between the two spheres; its value 
can be predicted analytically, and is given by the following 
equation: 

ptt/2 



On 



J* /Jl F(0)D( sirif9)df? 



J i ; /2 F(9)de 



where 



F(9) 



sin 3 8 cos# 



sin 2 0rf n /r, 



1/4 



(10) 



(11) 



and 
D(s) = 1 



a 2 



(1 - (s 2 /rLt)) - 04 (1 - (s 2 /r 2 out )) 2 . (12) 



The variation of the calculated value of the interior anguar 
velocity as a function of the width of the g ap — 7"out Tin 
is presented in Fig [j. It can be seen that the simulations, 
represented by the black squares, fit well the analytical pre- 
dictions provided that the Ekman number is small enough 
(i.e. below 10 -6 ). It is also interesting to note, as an aside, 
that for gap width of about 3% of the radiative zone's radius 
(which corresponds to the width of the solar tachocline), 
the interior angular velocity is 93% of the equatorial veloc- 
ity, which is very close to the value observed. It is not clear 
whether this interesting match is a mere fortuitous coinci- 
dence or the result of some more subtle physical processes. 

Another way of comparing the results of the simula- 
tions to analytical predictions is through the construction 
of the Ekman spiral, which is a parametric representation of 
the azimuthal velocity against the latitudinal velocity as a 
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Figure 2. Rotation profile and streamlines in the non-magnetic case for E v = 8 X 10 . The streamlines are line-coded so that clockwise 
motions are represented by dotted lines and anti-clockwise motions are represented by solid lines. The interior angular velocity in this 
simulation is Q ln = 0.75f2 C q- 
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Figure 4. Predicted and calculated Ekman spirals. The analyt- 
ical prediction in shown as a continuous line, whereas the sim- 
ulations are shown as squares for E v = 10 — 4 or triangles for 
E v = 6.5 X 10 — 6 . The inset shows an enlargement of the central 



function of radius, at a fixed co-latitude 6. Fig. ^| compares 
the results of the asymptotic solution and the true numeri- 
cal solution for a slightly different simulation, in which the 
angular-velocity profile imposed on the outer boundary is 
chosen to be constant with value fi out = Q c + 10~ 5 , and the 
angular velocity of the inner core is simply Q- m = f2 c (the 
no-torque condition is dropped). Note how the fit of the 
asymptotic analytical prediction to the numerical solution 
is valid only provided the Ekman number is small enough. 
For larger Ekman numbers, viscosity plays a non-negligible 
role in the dynamics of the fluid outside the boundary lay- 
ers, invalidating Proudman's asymptotic analysis. This re- 
sult has been obtained already by Dormy, Cardin & Jault 



(1998) who studied the non- magnetic case extensively. The 
good agreement between the analytical asymptotic solutions 
and the simulations validates the numerical procedure. 



4 RESULTS IN THE MAGNETIC CASE 

The influence of the magnetic field on the fluid depends es- 
sentially on two parameters: the field strength and the mag- 
netic diffusivity. In this section, three regimes are presented 
for varying Elsasser number at fixed (E v , E v ). The Elsasser 

the dependence of 



number is then fixed, and in Section 4.2 



the solution on the magnetic Ekman number is presented. 



4.1 Varying the field strength 

In these first simulations, only the Elsasser number A is 
varied. Note that the definition of A defined in equation 
(^) uses the value of the amplitude of the radial component 
of the magnetic field on the inner boundary. Accordingly, 
it should normally be defined using the true value of the 
density at r = ri n , which in reality is of order of pi n = 20 
g cm" 3 rather than the chosen uniform value of po = 1 g 
cm" 3 . As a result of the uniform-density approximation, the 
values of A chosen in this work should be interpreted with 
care and regarded as rough indications rather than precise 
values. Note also that as the "initial" magnetic field varies 
as 1/r 3 , the typical local Elsasser number near the surface 
is about 500 times lower than A. 

The values of the viscous and magnetic Ekman numbers 
chosen for these simulations are identical: E v — E v — 2.5 x 
10" 4 . This value was chosen for convinience, as it is then 
easy to obtain solutions for any value of the magnetic field 
strength. 



4-1.1 Low-field case, A = 1/25 

This first simulation is shown in Fig. |^, which presents the 
result in the case of a low Elsasser number (A = 1/25). This 
corresponds to Bo = 0.25 T. The local Elsasser number near 
the surface is of order of 7 x 10 . 
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(a) Rotation rate Q/Q 




(b) Streamlines 
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Figure 5. Simulation results for A = 1/25, E v = 2.5 X 10 -4 , and E n = 2.5 X 10~ 4 . Same line-style coding as in Fig. ti for the streamlines. 



The structure of the interior angular velocity is domi- 
nated by Coriolis forces, and the angular- velocity profile is 
close to Proudman (cylindrical) rotation, except maybe near 
the inner core where the influence of a magnetic field can 
be seen through the slight deviation in the angular-velocity 
contour lines. Because of the additional Lorentz forces in 
the momentum equation, the circulation is no more limited 
to cylindrical surfaces and takes a rather different pattern, 
with two cells that burrow deeply into the radiative zone. 
As expected from standard Ekman-Hartmann pumping (see 
Acheson & Hide, 1973), the typical value of the latitudinal 
component of the flow near the top boundary is of order 
of rQ., i.e. ug ~ 1.5 x 10 -2 in units of r c Q c . The advection 
of the poloidal field by the circulation can be seen in Fig. 
^|(d) where the field lines are dragged equatorwards in the 
interior and polewards near the top boundary by a strong 
anticlockwise polar cell. The shear persists throughout the 
fluid region, and as a result, leads to the winding up of the 
poloidal field into a relatively strong toroidal field. Typical 
values of the toroidal field are of order of one tenth of the 
value of the poloidal field near the core. This structure shows 
little resemblance with the observations, failing in particular 
to impose uniform rotation within the fluid region. 



4.1.2 High-field case, A = 25 

The second simulation is shown in Fig. [], which presents 
the result in the case of a high Elsasser number. This corre- 
sponds to Bo = 6.5 T. The local Elsasser number near the 
surface is typically of order of 4.5 x 10 -2 . 



In the strong-field case, (see Fig. g), the system is es- 
sentially dominated by the Lorentz forces everywhere but 
in diffusive boundary layers. The magnetic field is hardly 
affected by the circulation, and keeps essentially its dipo- 
lar structure everywhere within the fluid region. As a result 
of Ferraro's isorotation law, the constant angular-velocity 
contours follow closely the magnetic field lines. Magnetic 
connection with the inner core is strong, and a large region 
within the outermost closed field line near is forced to rotate 
nearly uniformly with angular velocity fl[ n . The matching 
of this uniformly rotating region to the polar regions and to 
the outer boundary occurs through a diffusive internal shear 
layer of thickness 5ii 

1/4 

<h = ( 1 , (13) 



Aloe 

where Ai oc is a local Elsasser number, defined with the local 
value of the typical amplitude of the field instead of Bo- A 
mo re detailed study of this shear layer is outlined in Section 
and presented by Kleeorin et al. (1997), Dormy, Cardin 



5.2 



& Jault (1998) and Dormy, Cardin & Soward (2001). 

In the polar regions, the field lines are anchored into 
the differentially rotating convection zone and the latitudi- 
nal shear is transmitted to the radiative zone through Fer- 
raro isorotation. The matching of the rotation profile in the 
polar regions with the outer boundary and the inner core 
is ensured by the existence of a diffusive Ekman-Hartmann 
layer with thickness 5± such that 

1/2 

5x = ( 1 • (14) 



Aloe 
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(a) Rotation rate Q/fi 



(b) Streamlines 
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Figure 6. Simulation results for A = 25, E v = 2.5 X 10 , and = 2.5 X 10 . Same line-style coding as in Fig. ti for the streamlines. 



This boundary layer is discussed in detail in Section j5.l| . 

The circulation is essentially limited to equatorial re- 
gions, with one strong principal cell following the internal 
shear layer, and weak secondary ones. Again, typical flow ve- 
locities along the shear layer are of the order of rQ, whereas 
velocities across the layer are of the order of 5iif2. This flow 
is however too slow to have any significant effect on the field. 
In that case, linear asymptotic analysis can be performed on 
the system, as shown by Kleeorin et al. (1997). The toroidal 
field is mostly limited to regions of shear (near the poles) 
with a very small amplitude in the co-rotating regions. It 
is worth mentioning that in this case, because the inner re- 
gions rotate almost uniformly, relaxing the rigidity condition 
within r = n n is likely to have little effect on the solution. 



4-1.3 Intermediate-field case, A = 1 

This third simulation is shown in Fig. ^, which presents the 
result in the case of an intermediate value of the Elsasser 
number. This corresponds to Bo = 1.3 T. The local Elsasser 
number near the surface is typically of order of 2 x 10 -3 . 

The intermediate-field case reveals the emergence of two 
distinct regions: in the interior the system is dominated by 
the magnetic field, and is in a state close to uniform ro- 
tation, with a large region rotating with angular velocity 
flin (except in a narrow latitude band around the polar re- 
gions). Most of the shear is confined to a thin layer, the 
"tachocline" ; however, within that shear layer and especially 
near the equator, the system can be observed to follow cylin- 



drical rotation, showing that the system is dominated by 
Coriolis forces rather than Lorentz forces. 

The following phenomenon is happening. A two-cell 
circulation with upwelling in mid-latitudes is driven by 
Ekman-Hartmann pumping near the outer boundary, as be- 
fore. However, in this case the typical circulation velocity 
(tie ~ rSl) is sufficiently high compared to the Lorentz 
stresses to advect the magnetic field in a direction parallel to 
the outer boundary in most latitudes, except near the equa- 
tor, where it is advected downwards, and near the poles, 
where the flow is parallel to the field. As a result of this 
advection process, the amplitude of the radial field on the 
boundary is everywhere reduced, which confines the field 
to the radiative zone and diminishes the magnetic connec- 
tion between the interior field and the imposed latitudinal 
shear. This can be observed particularly well in Fig. 0, which 
shows that the radial component of the field is close to zero 
in the "tachocline" region. Moreover the downwards advec- 
tion of the field in the equatorial regions reduces the value 
of the local Elsasser number (proportional to the square of 
the amplitude of the field), which explains the emergence 
of the cylindrical rotation in that area. Flux conservation 
implies that the magnetic field strength is correspondingly 
increased in regions just below. This magnetic field evacu- 
ation by the circulation can also be seen in Fig. |l(] which 
shows the square of the amplitude of the magnetic field on 
the equator as a function of radius. The plot shows particu- 
larly well the two regions: 

(i) in the core, the field is hardly perturbed by the differ- 
ential rotation imposed on the top, and varies with r -3 just 
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Figure 7. Radial component of the magnetic field as a function 
of the normalized radius at cola titude 9 = 7r/3, for the simula- 
tion presented in Section 4.1,3| (i.e. with the parameter values 



E v = 2.5 X II 
in Section (4.2 



E v = 2.5 X 10~ 4 and A = 1) (dotted line), and 
(i.e. with the parameter values E v = 6.25 X 10 -5 , 
E v = 6.25 X 10~ 5 and A = 1) (solid line). The dashed line rep- 
resents the same quantity for a non-rotating system, where the 
magnetic field solution is a dipolar field decaying as r -3 . 



Figure 10. Normalized square of the amplitude of the magnetic 
field on the equator as a func tion o f normalized radius for the sim- 
ulation presented in Section 1.1. s| (i.e. with the parameter values 



E„ = 2.5 x 10 
and Section 



2.5 X 10~ 4 and A = 1) (dotted line), 
4.2 (i.e. with the parameter values E v = 6.25 X 10~ 5 , 
E v = 6.25 X 10~ 5 and A = 1) (solid line). The dashed line rep- 
resents the same quantity for a non-rotating system, where the 
magnetic field solution is a dipolar field decaying as r~ . 



as the field would were the system not rotating (as repre- 
sented by the dashed line). 

(ii) the advection of the field by the circulation can easily 
be seen near the surface: just below the convection zone, 
the amplitude of the field is significantly smaller than the 
initial dipolar field, and slightly lower down the amplitude 
is correspondingly higher, as required by flux conservation. 

Note that the two regions are separated by a magnetic dif- 
fusion layer. Additional simulations for different viscous and 
magnetic Ekman numbers suggest that the thickness of the 
layer is, as expected, of order of 8\\ . 

Conversely, the magnetic field keeps the circulation 
from burrowing deep into the radiative zone and confines 
it to a shallow region. This confinement can be seen in Fig. 
^| but is represented best in Fig. [ll], which shows the latitu- 
dinal component of the velocity as a function of radius. Note 
how the circulation is heavily suppressed below r = 0.9r c . 
Finally, one can see that the polar regions, which are unaf- 
fected by the circulation, propagate the slow polar rotation 
at all radii into the radiative zone, as in the high-field case. 



4.2 Varying the magnetic diffusivity 

The amplitude of the magnetic field is now set to be that 
of the intermediate field strength, Bo = 1.3T. When the 
magnetic diffusivity is decreased by a factor of four, as it 
is shown in Fig. ^, advection of the magnetic field by the 
circulation becomes more important compared to diffusion. 
This has several consequences: 

(i) the fluid is closer to being a perfect fluid, driving the 



the system closer to isorotation. A larger volume of fluid 
in the interior is rotating nearly uniformly with the inte- 
rior angular velocity and the slow rotating polar regions are 
confined within a smaller latitude band; 

(ii) the field confinement and the reduction of magnetic 
stresses is more efficient, as can be seen in Fig. ^ and Fig. ^. 

(iii) as can be seen in Fig. ^ and Fig. ^ the magnetic 
evacuation in the surface equatorial regions is also greater, 
and occurs more abruptly as expected. The typical width 
of the transition layer between the magnetically dominated 
interior and the magnetic free tachocline seems to vary as 
(_E„_E I) ) 1//4 , confirming that it is a si mple magnetic diffusion 
layer of the type analysed in Section 



5.2 



(iv) as a result, the circulation is confined within a smaller 
volume by the magnetic field (see Fig. ^| and Fig. pd[ ). 

As mentioned previously, quantitative results of the 
model cannot provide any serious predictions of the obser- 
vations (because the inadequate driving mechanism used for 
the circulation, because of the assumption of incompressibil- 
ity and uniform density, and because of the gross discrep- 
ancy between the real solar diffusivities and the ones used 
in these simulations). Nonetheless, the calculated interior 
angular velocity provides an interesting point of compar- 
ison between different solutions corresponding to different 
parameters, and reflects on the various dynamical processes 
occuring in the radiative zone. Figure [13 shows the pre- 
dicted ratio of interior to equatorial angular velocities as a 
function of the Ekman number in the intermediate-field case 
and high-field case as a function of magnetic Ekman number. 
Again, the viscous Ekman number is chosen to have the same 
value as the magnetic Ekman number. The main reason for 
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(a) Rotation rate Q/fi B 



(b) Streamlines 




Figure 8. Simulation results for A = 1, E v = 2.5 X 10 4 , and E v = 2.5 X 10 4 . Same line-style coding as in Fig. La for the streamlines. 
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Figure 9. Simulation results for A = 1, E v = 6.25 X 10 5 , and E v = 6.25 X 10 5 . Same line-style coding as in Fig. bj for the streamlines. 
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the inner core, with a stronger connection for lower Ekman 
number. However, as the Ekman number is decreased past 
a certain threshold, shear expulsion from the core slowly 
starts affecting the polar regions and the inner core velocity 
increases again towards the observed values. Note that the 
value of the interior angular velocity in the extremely diffu- 
sive case is of 96% of the equatorial angular velocity, a value 
which was predicted by Gough (1985). 



5 COMPARISON WITH ASYMPTOTIC 
ANALYTICAL ANALYSIS 

This section focuses on presenting two asymptotic solutions 
near the spherical boundaries, in the case where the mag- 
netic field is mostly perpendicular to the boundary (which 
occurs near the poles on the inner boundary) and in the case 
where the magnetic field is mostly parallel to the boundary 
(which occurs on the outer boundary near the equator). The 
boundary layer solutions are then compared with the numer- 
ical solutions. 



Figure 11. Latitudinal component of the velocity as a function 
of the normalized radius at colatitude 8 = n/3, in units of the 
azimu thal velocity r c Sl c , for the simulation presented in Section 



4.1.3 



(i.e. with the parameter values E v 



2.5 X 



2.5 X 10 4 and A = 1) (dotted line), and in Section 4.2 (i.e. with 



la 



the parameter values E v 
A = 1) (solid line). 
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Figure 12. Interior angular velocity as a function of the Ekman 



number E - 



E v for the low field case (triangles), interme- 



diate field case (stars) and the high field case (squares) 



the initial decrease in Q ln with decreasing Ekman number 
is the following: for this range of magnetic Ekman numbers, 
the polar regions are hardly affected by the shear expulsion, 
and rotate with a low angular velocity more or less at all 
radii. A slowly rotating region, connected to the core via 
magnetic stresses, necessarily imposes its slow rotation to 



5.1 Analysis of the boundary layer in the polar 
regions 

5.1.1 Analytical derivation of the boundary layer solution 

When the magnetic field is mostly perpendicular to the 
boundary, an Ekman-Hartmann boundary layer develops 
(see the review by Acheson & Hide (1973)). A derivation 
of the boundary layer analysis is now repeated for the sake 
of clarity. Assuming that the magnetic field is essentially 
perpendicular to the boundary with constant amplitude Bo, 
one can write 



B = B e r 



(15) 



Calling x = r — rjn, the simplified system of MHD equations 
in the boundary layer on the inner sphere is 

^dS ^ d 2 L 
A— — h E v —— 

ox ox- 1 



9 dip 



i_&p 

rin dx 
9 dL 



db 

. d 2 b ^ d A i> 



dL 

dx 



-7T- = ~E, 



d 2 s 



(16) 



where ip is the stream function of the fluid flow, such that 



smdug = - — 
r or 



r 2 dfi 



and 



L = r sm8u^, , S = r sinOB^ and b = sin^Bg . 
Grouping these equations yields 



F 



(17) 



(18) 



(19) 



where F represents 
lution of the kind F oc e 1±x yields 



d 2 S d 2 b d 2 f 
dr 2 



or ^p. Looking for a so- 



- = (A - E u E v yl 



2\2 



(20) 
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and in turn, that fj_ = ±(/3 ± ia) where 

1/2 



A + ^A 2 + 4^2^ 



2i?„ E v 



A 



E V E V 



1/2 



-A+yA 2 +V£g 



1/2 



«/3 



(21) 



for small enough _E„ and E 1 ,,. The solutions, which must be 
bounded, can then be written out as 



Q ~ Q c e"' 33: + Q 



(22) 



where Q is either of the quantities S, b, L or ip, and Q c 
and Qo are integration constants. Note that the boundary 
layer thickness 8± — 1//3 is, as expected, the same as the 
one obtained by MacGregor & Charbonneau (1999) in their 
open-field calculations. 



0.005 



-0.01 



-0.015 



5.1.2 Comparison with the numerical solutions 

In order to compare rigorously the simulations to the ana- 
lytical solutions derived above, the matching of the solutions 
obtained in the boundary layer to those in the bulk of the 
fluid should be performed. However, the solution in the bulk 
of the fluid, in particular in the polar regions, is dominated 
by geometric effects (as the latitudinal derivatives are not 
necessarily negligible) and diffusive effects (as the Ekman 
numbers used in the simulations are not small enough to 
justify neglecting the diffusive terms); as a result, it is be- 
yond to scope of this analysis to derive an analytical solution 
for the solutions in the bulk of the fluid, and therefore at- 
tempt such a matching. However, it is still possible to check 
the qualitative behaviour of the solutions in the boundary 
layer. Fig. ^ shows the angular momentum function L as 
a function of the scaled variable £ = f3x — /3(r — n n ), for 
three different values of the Ekman numbers. Although the 
boundary conditions and asymptotic conditions are different 
for each case, this plots illustrates that the angular momen- 
tum function L indeed behaves as e l3x within the boundary 
layer (i.e. for values of £ below unity). 

5.2 Analysis of the boundary layer in the 
equatorial region 

5.2.1 Analytical derivation of the boundary layer solution 

Near the equator, the asymptotic analysis described in Sec- 
tion [0] breaks down, as the magnetic field is advected by 
the circulation into a a direction parallel to the surface. An- 
other type of boundary layer appears, which is analysed in 
this section. 

Assuming little variation in the latitudinal direction, 
the magnetic field is approximated by B m B^eg; also, as 
only the region near the equator is considered, one can set 
/i ~ in the following analysis. In that case, and in the 
boundary layer only, the MHD equations can be simplified 
to 

B dS 



-0.02 



A 



r D ut d/j, 



E V L 



Bo dL ,i 
r ou t on 



(23) 



10 



Figure 13. Variation with the scaled variable £ of the angular 
momentum function L at fixed co-latitude 9 = w/12 and fixed 
Elsasser number A = 1, for three different values of the Ekman 
numbers: E v = E v = 1CT 3 (solid line), E v = E v = 2.5 X 1CT 4 
(dashed line) and E u = E v = 6.25 X 10~ 5 (dotted line). 



These two equations entirely determine the variation of L 
and S near the boundary, and can be combined into 



d 2 Q 
dfj? 



2 E V E V d Q 

r out ' 



ABq <9C 4 



(24) 



where Q is either L or S, and where the boundary layer co- 
ordinate £ = r out — r was introduced. On the boundary, the 
latitudinal variation of L is given by the boundary condi- 
tions, so that 



-2(o 2 + l)fto 



-2(a 2 + 1)L 



So finally, the governing boundary layer equation is 
„/ | -i \ r 2 E v E n d i L 

~ 2(a2 + 1)L = r °^TBjW ' 



(25) 



(26) 



and similarly for S. The solutions which are bounded as 

£ — > oo are 



L = L c (p)e^ ( cos( 7 ||C) +L a ( M )e^ c sin( 7 ||C) 



and similarly for S, with 



7|| = 



2(a 2 + l)r 2 out E v E v 



1/4 



A, 



1/4 



(27) 



(28) 



where the primes denote differentiation with respect to r. 



Note that in the work of Riidiger & Kitchatinov (1997), the 
poloidal field is entirely confined within the radiative zone 
and, as a result, is everywhere parallel to the surface near the 
boundary. This type of boundary analysis therefore applies 
also to their results, with the following modification: as their 
latitudinal field varies as (1 — C) q l near the boundary, it can 
be shown that the thickness of the magnetic diffusion layer 
scales as 711 oc (A loc /^^) 1/(2+29) .The case presented here 
is recovered with q = 1. 
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Figure 14. Variation with normalized radius of (L — L c )/ (L eq — 
L c ) on the equator. The solid lin e is deduced from the numeri- 



cal solution presented in Section 1.2 



and t 



le dotted line is the 



exponential solution described in equation (2! 



5.2.2 Comparison with the numerical solutions near the 
equator 

Fig. jTi] presents the variation of (L — Li n )/(L cq — Lj n ) on 
the equator as a function of normalized radius, both for the 
numerical solution and for the "analytical prediction" given 
by 



Lin 



,711 (' — r oiit) 



Lin 



(29) 



where Lj n = rf n fl ln . This simulation corresponds to the pa- 
rameters E v = E„ = 6.25 x 10~ 5 and A = 1. If higher Ekman 
numbers are chosen, the asymptotic analysis does not satis- 
fyingly reproduce the solution. One can see that the general 
shape of the solution is well represented by the exponentially 
decaying solution with boundary layer width 5m = I/711. 

The purpose of this section was to determine an asymp- 
totic solution to the equations in the boundary layers near 
the inner and outer boundaries. This analysis has been used 
to optimize the mesh-points allocation procedure, as well 
as to check the behaviour of the numerical solution in the 
boundaries. 



6 DISCUSSION OF THE RESULTS 

6.1 Summary of the results and comparison with 
previous work 

Section ^ studied the interaction of a large-scale magnetic 
field with fluid motions in the radiative zone when a shear is 
imposed through magnetic and viscous stresses by the con- 
vection zone. Three regimes were studied, which depend on 
the intensity of the magnetic field. It was shown, by varying 
the magnetic field strength for a given set of diffusion pa- 
rameters (E V ,E V ), that there exists only a narrow range of 



parameters A that leads to a supression of the shear in the 
bulk of the radiative zone. When the system is dominated 
by Coriolis forces (A <C 1), the shear extends well into the 
radiative zone through the Proudman constraint. Moreover, 
two large circulation cells are driven by Ekman-Hartmann 
pumping on the boundary, and largely mix the radiative 
zone. In the opposite situation (A ^> 1) the system is domi- 
nated by Lorentz forces, and the effects of meridional motion 
are negligible. In that case, the magnetic field remains es- 
sentially dipolar and the latitudinal shear imposed by the 
convection zone extends along the field lines into the radia- 
tive interior. This system was studied in detail by Kleeorin et 
al. (1997) and Dormy, Cardin & Jault, (1998), who showed 
the existence of an internal shear layer along the outermost 
closed field line "associated with the recirculation of electric 
currents induced in the Hartmann layer as the internal Stew- 
artson layer is associated with the recirculation of meridional 
flows generated within the Ekman layers" . The thickness of 
this internal layer is of order of 5m (Dormy, Jault & Soward, 
2001). Note that in this strong-field case it is clear that the 
internal rotation profile depends entirely on the structure of 
the imposed field. This can be deduced from the comparative 
work of MacGregor & Charbonneau (1999), who show how 
an open-field structure, where the field lines are anchored in 
the convection zone, results in a differentially rotating ra- 
diative zone, whereas in a confined-field structure the shear 
can be confined to a thin tachocline. 

Finally, for A ~ 1, it is observed that shear exclusion 
indeed takes place. In that case, meridional motions signifi- 
cantly distort the poloidal field structure, to confine it grad- 
ually to the radiative interior as it was initially suggested by 
Gough & Mclntyre (1998). This confinement results in the 
reduction of the magnetic stresses connecting the radiative 
interior to the differentially rotating convection zone, and 
to the expulsion of the shear by the interior field to a thin 
shear layer, as in the simulations first presented by Poidiger 
& Kitchatinov (1997). However, there exists a fundamen- 
tal difference between the shear layer (tachocline) thickness 
and position predicted from the Poidiger & Kitchatinov and 
MacGregor & Charbonneau (1999) simulations and the ones 
presented here. In their model for the confined field, as the 
field structure is prescribed to the system and the meridional 
flows are neglected, the position of the magnetic diffusion 
layer is at the boundary and its width is exactly propor- 
tional to 5ii(So), where Bo is the imposed value of the lati- 
tudinal field at the boundary. In the simulations presented 
here, it is shown that the meridional motions can affect sig- 
nificantly the strength and the geometry of the field near 
the top of the radiative zone through advection processes, 
which controls the position and thickness of the magnetic 
diffusion layer. This phenomenon is intrinsically linked with 
the nonlinear interaction of the field and the flow and can be 
studied only through the resolution of the nonlinear MHD 
equations. The influence of the flow on the field depends also 
on the strength and on the geometry of the flow, which in 
turn are intrinsically related to the mechanism which drives 
the flow. These comments naturally build up to the model 
proposed by Gough & Mclntyre (1998). 

One of the principal features of the Gough & Mclntyre 
model is their suggestion that the magnetic diffusion layer 
could be situated below a virtually magnetic-free tachocline 
(by magnetic-free, they refer more to a system where the in- 
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fluence of the magnetic field in the dynamical equilibrium is 
negligible than to a system that has no magnetic field at all) . 
This structure would result from the almost complete con- 
finement of the magnetic field within the radiative zone by 
the meridional motions (except perhaps in a localized up- 
welling region in mid-latitudes). Although the simulations 
above cannot reproduce quantitatively the dynamical bal- 
ance proposed by Gough & Mclntyre (as the circulation in 
this paper is artificially driven by Ekman-Hartmann pump- 
ing instead of thermal imbalance) , these qualitative features 
can easily be identified in Fig. [?] and Fig. |Io| . The system 
seems indeed to lead naturally to a magnetically dominated 
interior and a rotationally dominated tachocline separated 
by a magnetic diffusion layer of a thickness that is observed 
to vary as 5n . 



6.2 Comparison with observations 

The next step of this analysis consists in comparing the nu- 
merical simulations with the observations. In this compar- 
ison process, it is essential to keep in mind two essential 
discrepancies between the model and the real sun: the driv- 
ing mechanism for the meridional flow in these simulations 
is artificial, so that the typical flow velocities may be er- 
roneous, and the typical Ekman numbers of the simulations 
are several orders of magnitude larger than in the sun. In the 
region of the tachocline, assuming that the flow is not tur- 
bulent, the magnetic and viscous diffusion coefficients are 
of order of v = 10cm 2 s _1 and r\ — 2 x 10 3 cm 2 s~ 1 , which 
implies that 



E v ~ 10" b and E„ 



2 x 10" 



(30) 



The main consequence of these discrepancies is that al- 
though the principal features of the interaction between 
fluid motions and large-scale fields in the sun can be studied 
through this method, it cannot provide reliable quantitative 
estimates. 

Despite this obvious shortfall, some aspects of the ob- 
servations can be reproduced in the simulations. 

(i) For sufficiently low magnetic diffusivity, a large region 
of the radiative zone is forced to rotate uniformly with the 
angular velocity of the core; such a uniform angular- velocity 
profile in the radiative interior is indicated clearly by the ob- 
servations. The simulations suggest that polar regions seem 
to rotate with a low angular velocity down to a latitude of 
about 50° for E v = E v — 2.5 x 10~ 4 , and down to lati- 
tudes of about 60° for E v = E v = 6.25 x 10" 5 : the polar 
shear is gradually confined to higher and higher latitudes as 
the magnetic Ekman number is reduced. Helioseismic inver- 
sions still have too low a resolution to provide any reliable 
observations of the polar regions; this numerical model how- 
ever predicts more slowly rotating fluid in the polar regions 
deep within the radiative zone, which is a natural feature of 
the dipolar field configuration (Riidiger & Kitchatinov 1997, 
Gough & Mclntyre 1998, MacGregor & Charbonneau 1999). 
Since the most likely configuration for a fossil field is dipo- 
lar, this polar feature is a robust prediction of this family of 
models. 

(ii) A thin boundary layer appears in which most of the 
shear is contained, particularly in the equatorial regions. 
The thickness of this shear layer is much larger than the 



observed thickness of the tachocline, a discrepancy which is 
again simply related to the high diffusivities adopted here. It 
is interesting to note, however, that the position of the mag- 
netic diffusion layer predicted by the simulations is deeper 
in the equatorial regions than at higher latitudes (which is 
due to the downward advection of the field near the equa- 
tor); this could be related to observations of the prolateness 
of the tachocline (Charbonneau et al, 1999). 

(iii) A two-cell meridional circulation is driven below the 
convection zone by Ekman-Hartmann pumping; this circu- 
lation is confined by the large-scale poloidal field to the 
tachocline region. This result validates the analysis of the 
sound-speed profile performed by Elliott & Gough (1999) 
and also relates to the upper limits in the depth of the 
tachocline mixing suggested by observations of the abun- 
dances of light elements in the convection zone (Brun, Turck- 
Chieze & Zahn (1999)). Again, the simulated depth of pen- 
etration of the circulation is much larger than suggested by 
the observations, a discrepancy which is due to the large 
diffusivities used for the numerical analysis, and to the as- 
sumption of incompressibility: comparison between high- 
and low-diffusivity cases indeed show a significant reduc- 
tion of the mixed-layer depth for lower Ekman numbers, and 
when the background stratification is taken into account, it 
is observed to impede radial motions dramatically(see sub- 
sequent work). 



Finally, as mentioned in Section 1.2, a useful point of 
comparison between different theoretical models and be- 
tween models and the observations is the value of the in- 
terior angular velocity fii n . Gough (1985) showed that for 
an incompressible tachocline controlled by viscous effects, 
Qin = 0.96f2 cq . MacGregor & Charbonneau (1999) pre- 
sented simulations in a confined field configuration which 
show a virtually uniform rotation profile for the radiative 
zone with fl[ n — 0.97f2 oq . The discrepancy with the observed 
value (fi c = 0.93f2 eq) is significant, and comparison with the 
Gough & Mclntyre model shows that either meridional mo- 
tions and the effects of heat transport and stratification are 
dominant in the dynamics of the tachocline, or that such a 
class of models cannot reproduce the observations and that 
very different dynamics are in play (as in the turbulent stress 
model proposed by Spiegel & Zahn, 1992, or the nonlinear 
development of MHD instabilities presented by Cally, 2001 
The prediction for the angular velocity presented in Fig. h 
are not conclusive when comparing them with observations 
as any discrepancy could be attributed to the large diffusiv- 
ities used in the simulations. 



6.3 Discussion of the model 

The ultimate aim of this work is to develop a self-consistent 
dynamical model of the tachocline which can be used to 
explain the large range of observations available. In order 
to do so, the idea proposed by Gough & Mclntyre is gradu- 
ally implemented into a numerically solvable nonlinear MHD 
model. Although previous work on those lines (Riidiger & 
Kitchatinov, 1997, MacGregor & Charbonneau, 1999) has 
already investigated some of the aspects of the interaction 
of a large-scale field and differential rotation, two essential 
ingredients to the model remained to be studied carefully: 
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(i) the nonlinear interaction of a meridional circulation 
and the imposed magnetic field 

(ii) the baroclinic driving of meridional motions and the 
effects of compressibility and stratification. 

This first paper focuses on the first point only by artificially 
driving a meridional flow with Ekman-Hartmann pumping 
on the boundaries. It is found that this artificial system re- 
produces qualitatively (but not quantitatively) most of the 
aspects of the Gough & Mclntyre model, and of the obser- 
vations. The good qualitative agreement with the Gough & 
Mclntyre model is perhaps surprising, but is simply due to 
the qualitative similarities between the Ekman-Hartmann 
flow and the baroclinicly driven flow, in particular with 
respect to the two-cell structure with upwelling in mid- 
latitudes. In both cases, the quadrupolar structure of the 
flow is linked with the transition in mid-latitudes from a neg- 
ative radial shear to a positive radial shear. When Ekman- 
Hartmann layers are present on the outer boundary, they 
naturally drive poleward motions at high latitudes and equa- 
torward motions at low latitudes, with a return flow up- 
welling in mid-latitudes. In the sun, thermal- wind balance 
across the tachocline is shown to lead to temperature fluctu- 
ations with a minimum in mid-latitudes, which also results 
in an upwelling in this region (Spiegel & Zahn, 1992). 

Another important point concerns one of the principal 
assumptions of this numerical model: the linearization of the 
problem with respect to the advection of momentum by the 
meridional flow. In order to verify th at ( u ■ V)u can indeed 
be neglected, it is compared in Fig. ha to the linear term 
2e z x u as well as th e Lo rentz force j x B, for the simulation 
presented in Section 4.2 and in Fig. ^. As assumed, the non- 
linear advection term is everywhere much smaller than the 
Coriolis term, except perhaps in the Ekman-Hartmann lay- 
ers near r — n n where both Coriolis and advection term van- 
ish (the system is then dominated by the magneto-viscous 
balance). It is therefore justified to neglect it. A similar com- 
parison can be performed (see Fig. |l5|) between the vorticity 
advection term Vx(u-Vu), the background vorticity advec- 
tion term 2Vx (e z x u) and the Lorentz stresses Vx(j'x B); 
the conclusion is the same. 

Finally, a few more points must also be raised. 

(i) As mentioned previously, the diffusivities used in the 
simulations are far higher than in the sun. The principal ef- 
fect of high diffusivities is to reduce the interconnection be- 
tween the field and the flow which can significantly change 
the balance of forces within the fluid region. This discrep- 
ancy is forced by numerical solvability and can be reduced 
at great cpu expense only. 

(ii) The boundary conditions on the magnetic field are 
found to have significant influence on the problem (MacGre- 
gor & Charbonneau 1999). Although the systematic study 
of the effects of boundary conditions on the system has not 
been discussed in this paper, it is illustrated, for example, in 
the comparison between this work and the results of Dormy, 
Cardin & Jault (1998) who study (in a geophysical applica- 
tion) the effects of an insulating outer boundary; as a result, 
the structure of the toroidal field in the interior is altered 
significantly. 

(iii) The boundary conditions near the top boundary 
(which is the base of the convection zone) poorly repre- 
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Figure 15. Lower panel: comparison between the advection term, 
the Coriolis force and the Lorentz force at co-latitude 9 = 60° 
in units of r c Q^, for the simulations presented in Section 4.2. 
Upper panel: comparison between the vorticity advection term, 
the background vorticity advection and the magnetic stresses at 
co-latitude 9 = 60° in units of f2;?. 



sent the dynamical interaction between the flow and field 
in the tachocline and the convection zone: the convection 
zone, and to some extent the overshoot region below, is the 
seat of fully developed turbulence. It does not behave as an 
impermeable wall to the circulation in the tachocline, or as a 
uniformly conducting medium to the magnetic field. More- 
over, it is for example widely believed that the interface of 
the tachocline and the convection zone is the seat of the solar 
dynamo, which is a non-steady magnetic structure oscillat- 
ing with a 22 yr period. How can this structure be matched 
with the assumed steady, nearly dipolar field of the interior? 
The problem was avoided in this work by assuming that the 
upper boundary of the tachocline lies safely below the re- 
gion of influence of the dynamo (see Garaud, 1999) and the 
overshoot region, but this is clearly an over-simplification. 

(iv) This also raises one of the major issues of the dy- 
namics of the tachocline: to what extent can a steady model 
represent the system? Firstly, how can the interior field be 
represented by a steady field when part of it diffuses, or is 
advected away into the convection zone? In the Gough & 
Mclntyre model, the diffusion of the magnetic field out of 
the radiative zone is impeded by the confining action of the 
circulation on the poloidal field, which supports the steady- 
state assumption. However, the upwelling region discussed 
by Gough & Mclntyre (and observed in Fig. ^) drags the 
poloidal magnetic field into the convection zone, where it 
could occasionally reconnect. This phenomenon is intrin- 
sically non-steady, although possibly quasi-periodic (as no 
magnetic flux is lost through reconnection) . Secondly, the 
hydrodynamical stability of the solution to various insta- 
bilities (and in particular to non-axisymmetric perturba- 
tions) should be assessed. It was shown by Garaud (2001a) 
that two-dimensional hydrodynamical instabilities are un- 
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likely to be strong enough to have a significant effect on 
angular-momentum redistribution in the tachocline. How- 
ever, the stability of the shear flow tends to be upset by 
the presence of a magnetic field, as suggested by Gilman & 
Fox (f997), Dikpati & Gilman (f999) and more recently by 
Cally (2001). In a stratified fluid, the buoyancy of toroidal 
flux tubes may lead to a magneto-convective instability. Fi- 
nally, Alfven waves can propagate along the magnetic field 
and lead to the presence of either "localized" oscillations, 
or global torsional oscillations. The combined effect of these 
instabilities and waves on the flow may be important in re- 
distributing angular momentum within the radiative zone 
and across the tachocline. This has not been taken into ac- 
count in this model. 

In future work each of these problems must be ad- 
dressed. The gradual convergence of the system towards 
lower and lower diffusivities is "merely" a computational 
challenge, which can in principle be solved. There is hope 
that for low enough diffusivities the system reaches an 
asymptotic state which is more-or-less independent of the 
Ekman parameters. Asymptotic analyses, following the steps 
of Kleeorin et al. (1998) for instance, may provide a way of 
by-passing boundary layers in the flow by providing jump- 
conditions across the layers. This would greatly help reduc- 
ing some of the numerical difficulties. However, looking at 
the importance of geometric or nonlinear effects in the sys- 
tem, it is unlikely that an asymptotic analysis could provide 
a full solution of the problem. The problem of the represen- 
tation of the dynamical interaction of the tachocline and the 
convection zone through adequate boundary conditions is a 
strong theoretical challenge, which will be looked into in the 
future. The most obvious route is through the inclusion of 
the convection zone in the computational domain, and the 
ad-hoc prescription of Reynolds stresses which would lead 
to the observed rotation profile. But comparing the relative 
importance of all the points mentionned in this section, it 
is clear that the next step towards improving the model is 
through the inclusion of the effects of compressibility, energy 
transport and stratification on the dynamics of the system, 
which will allow a correct quantitative representation of the 
baroclinic driving of the meridional motions as well as the 
effects of the strong stratification in the radiative zone. This 
is work which will be presented in a subsequent paper. 



7 CONCLUSION 

This paper presents a numerical analysis of the nonlinear 
interaction between a primordial field and large-scale fluid 
flows in the solar radiative zone when a latitudinal shear is 
imposed from the convection zone through a combination of 
magnetic and viscous stresses. 

Within the scope of some simplifying assumptions (ax- 
isymetry, steady-state, incompressibility), the presence of a 
large-scale field in the radiative zone was shown to lead to 
three dynamical regimes, which depend essentially on the 
strength of the underlying field. When the Elsasser number 
is low, the flow is dominated by the Proudman constraint, 
which enforces a rotation profile roughly constant on cylin- 
drical surfaces throughout most of the radiative zone. When 
the Elsasser number is high, the magnetic field propagates 



the shear to a large part of the interior according to the Fer- 
raro isorotation law, and only within the outermost closed 
field line does one observe uniform rotation. When the El- 
sasser number is of order of unity, the magnetic field is weak 
enough to be advected by meridional motions driven by 
Ekman-Hartmann pumping on the boundary. This advec- 
tion process stretches the field along the boundary (thereby 
reducing the magnetic stresses connecting the convection 
zone and the radiative zone) and pushes the poloidal field 
deeper into the interior. As a result, the shear is confined to 
some shallow layer below the convection zone, which can be 
likened to the tachocline. Conversely, the meridional circu- 
lation is kept from burrowing into the radiative zone by the 
accumulation of magnetic flux below the tachocline. 

The structure of the solution in the intermediate-field 
case is in qualitative agreement with the angular- velocity ob- 
servations, as well as the predictions of the Gough & Mcln- 
tyre (1998) model. The simulations also show the presence of 
efficient mixing in the tachocline which would influence the 
abundances of chemical elements in that region; this can also 
be detected observationally (Elliott & Gough, 1999, Brun, 
Turck-Chieze & Zahn, 1999). 

As expected, the quantitative agreement between the 
model and the observations is poor; this is due on the one 
hand to the large diffusivities used in the simulations, and on 
the other hand to some oversimplification of the dynamics 
involved (through the boundary conditions chosen, and also 
the assumption of incompressibility, which drive a merid- 
ional flow quantitatively different from what one could ex- 
pect in the tachcoline). Both of these problems will be ad- 
dressed in future work; in particular, the effects of compress- 
ibility will be studied in a following paper. 

In conclusion, although the dynamical behaviour of the 
sun is likely to be far more complex than that of the sim- 
ulations presented in this paper, it is clear that these sim- 
ulations can grasp, and explain, some of the fundamental 
aspects of the dynamics of the observed solar differential 
rotation. Moreover, they lay a strong basis for future im- 
provements of the models along the lines described above. 
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APPENDIX: ASYMPTOTIC SOLUTION IN THE 
NON-MAGNETIC CASE 

Following the work of Proudman (1956), the equations 



zeroth order in E v 



1/2 



[e z x u)^ = E„(V 2 u) 4 
[Vx(e z x u)L = E v {V 2 u) 4 



(31) 



are first solved in the main body of the fluid, then succes- 
sively near the bottom and top boundaries. In the main body 
of the fluid, the viscous stresses are negligible, so that u is 
independent of z, and the solutions are 



i/> = tpo(s) and X = Xo(s) 
where \ i s defined such that 



X 



r smB 



(32) 



(33) 



In order to solve the problem near the lower boundary, 
Proudman introduces the stretched variable £ such that 



C = (r - r in )E- 1/2 cos 1/2 6 



(34) 



This recognizes the presence of a boundary layer with thick- 
ness b~ v — eV 2 cos~ 1 ^ 2 9. The equations (pl|) become, to 



= -4 



dip 



9C 
dx 



(35) 



<9C 5 

dhc 

<9C 5 "dC 

Note that these approximations are not valid near the poles, 
where latitudinal derivatives may become important, and 
near the equator, where the thickness of the boundary layer 
S diverges. 

Assuming that the whole system is rotating with angu- 
lar velocity f2i n so that effectively fi c = ^in (this only calls 
for a slight re-definition of the Ekman number), matching 
with the boundary conditions at the bottom boundary re- 
quires that x ~ * as r — ► 1, or £ — » 0. Moreover, the im- 
permeable boundary condition requires that ip — on the 
boundary, as well as dip/d^ = 0. The solution to equations 
( pB[ ) which fulfills all these boundary conditions, and which 
is bounded as £ — > oo is 

V(C,0) = ^i(e)(l-e- c (cosC+ sinC)) 

X(t,0) = Xl(6) (l-e~ c cosC) , (36) 

where tpi(9) remains to be determined, and 

X i{0) = 2E~ 1/2 cos 1/2 0V>i(0) • (37) 

As £ — > oo, these functions must match onto the solution 
obtained previously for the main body of the fluid region so 
that is then easy to see that one must have 

Xi(0) = Xo(nn sm9) and V>i(0) = M^n sm6) . (38) 

This result can be combined with equation (^) and yields 
the matching condition 

Xo(r ln sin0) = 2£ , ~ 1/2 cos 1/2 (97/>o(r in sin(9) . (39) 

In order to study the boundary layer near the top boundary, 
another stretched variable is introduced: 



)E' 1 ' 2 



Mil 



(40) 



The scaled equations are the same as before (cf equations 
(pB|)); the boundary conditions for the stream function arc 
also the same as for the lower boundary when x ~ * 0, but 
the differential rotation must now match onto that of the 
convection zone, so that 



Xi r = fout,0) = r ont sin 6Q CZ (0) 



where 



Clcz{9) = fi cq (l — 0,2 COS 2 9 — £14 cos 4 (9) — Qi n 



(41) 



(42) 



The solutions to equations 



which fulfill these conditions 



= (l-e^( cos£+ sinO) , 

X(£,0) = X2(9) + 2E- 1/2 cos 1/2 ^ 2 (6»)e-« cos^ , (43) 
with 

X2(9) = r 2 ut sin 2 9Q cz (9) - 2E~ 1/2 cos 1/2 '0^,(0) . (44) 

As before, matching with the solution in the main body of 
the fluid implies that 

i>2(9) = ipo(r ut sinS) and X2{9) = Xo^out sin$) , (45) 
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so that 



Xo(r ou t sin6») = r out sin 0Q CZ (0) - 2_E„ cos Oipo(r ou t sin6 

(46) 

Since i/)o an d Xo are functions of s only, the two match- 
ing conditions given by equations (^) and ([Ifi]) can also be 
rewritten as 

Xo (s)^2E- 1/2 (l-(s/r in ) 2 ) 1/4 Ms) , 
X0 (s) = » 9 J4(a) - 2£" 1/2 (1 - ( S /r out ) 2 ) 1/4 ^o( S ) , (47) 
where 

£l' cz {s) = fJoq [l - «2 (l - (s/r out ) 2 ) 

-04 (1 - (s/r out ) 2 ) 2 ] -n in , (48) 
which can now be solved uniquely as 

El' 2 s 2 nUs) 



ipo(s) =- 



2 (1 - (s/r in ) 2 ) 1/4 + (1 - (s/r ou t) 2 ) 1/4 



Xo(s) * 2 (i-(*Mn) 2 ) 1/4 fU*) 

X ° [ ) (l-( S /n n ) 2 ) 1/4 + (l-( S Aout) 2 ) 1/4 ' 1 J 

The flow within two spheres is now known analytically ev- 
erywhere. 

Assuming that the sun is in equilibrium, the total 
torque applied by the convection zone on the radiative in- 
terior should be equal to that exerted by the solar wind on 
the convection zone. Since that torque is extremely small, 
it is assumed to be null as a first approximation, which is 
equivalent to requiring that the sun be in a steady state. 
This condition determines the interior angular velocity f2i n 
uniquely. 

In the non-magnetic case, the torques applied by the 
tachocline onto the radiative zone are purely viscous. As a 
result, the steady-state condition can be rewritten as 



T v {r = 1) = 2-kv 



3 . 2 a dQ 
r sm a— — 

or 



sin6»d6» = , (50) 



where T„ is the total viscous torque and £1 — fitn + 
(xA 2 sin 2 (9) is the total angular velocity at the base of the 
fluid region. Using the results derived previously, this con- 
dition can be rewritten as: 



n in _ j; /2 F(6)D( smd)de 



o= a r' 2 



F(d)A6 



(51) 



where 



= ^Llcose 

cos 1 /2e+(l- sin 2 6>r 2 n /r 2 ut ) 7 

and 

D(s) = 1 -oa (1 - (s 2 /r 2 ut )) -*M (1 - (s 2 /r 2 ut )) 2 . (53) 
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